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Abstract 

Public  domain  computer  programs  were  used  to  attempt  an  improved  model  of  the 
tritium  plume  observed  during  Macrodispersion  Experiment  2  (MADE- 2 ) ,  a  field 
scale  natural  gradient  experiment  conducted  at  Columbus  Air  Force  Base, 
Mississippi.  The  program  Geo-EAS  used  head  and  hydraulic  conductivity  data  at 
a  relatively  small  number  of  irregularly  spaced  test  locations  to  estimate 
corresponding  values  at  the  more  numerous  nodes  of  a  computational  grid  having 
66  rows ,  21  columns,  and  9  layers.  The  finite  difference  program  MODFLOW  was 
used  to  simulate  the  flow  of  groundwater  through  a  330  m  x  105  m  computational 
domain.  The  recent  BCF2  subroutine  package,  which  permits  rewetting  of  cells, 
allowed  the  vertical  discretization  to  be  more  accurate  than  in  previous  studies. 
Solutions  for  the  468  day  experiment  were  obtained  using  a  Sun  Sparcstation  2  for 
several  choices  of  convergence  and  storage  parameters.  The  simulations  had  small 
mass  balance  errors  and  were  consistent  with  continuous  head  observations.  The 
smallest  storage  coefficients  gave  the  best  agreement.  One  persistent  feature 
of  the  predicted  head  field  was  a  tendency  for  the  head  to  decline  -toward  the 
northwest.  This  suggests  that  the  plume  should  bend  toward  the  northwest,  but 
the  observations  show  a  bend  toward  the  northeast.  This  discrepancy  is  probably 
due  to  inaccurate  head  boundary  conditions  resulting  from  a  lack  of  piezometers 
in  the  northern  part  of  the  computational  domain.  The  flow  model  is  about  as 
accurate  as  the  data  permit. 

Tritium  plume  simulations  used  the  mixed  Lagrangian- Euler ian  finite  difference 
program  MT3D  to  solve  the  contaminant  transport  equation  using  the  MODFLOW- 
predicted  flow  field.  Thirteen  runs  were  made  using  various  advection  algorithms 
and  dispersivities ,  but  none  was  successful.  Numerical  instabilities  or  grossly 
unrealistic  predictions  ended  every  run  by  simulation  day  141.  Further  work  is 
needed  to  obtain  a  satisfactory  plume  prediction. 
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IMPROVED  NUMERICAL  MODELING  OF  GROUNDWATER  FLOW 
AND  TRANSPORT  AT  THE  MADE -2  SITE 

Donald  D.  Gray 

Dale  F.  Rucker 


INTRODUCTION 


Faced  with  the  need  to  remediate  groundwater  pollution  at  many  of  its  bases,  the 
Air  Force  has  undertaken  an  extensive  program  of  research  on  subsurface 
contaminant  transport.  The  Macrodispersion  Experiment  2  (MADE-2) ,  conducted 
together  with  the  Electric  Power  Research  Institute  and  the  Tennessee  Valley 
Authority,  was  a  key  element  of  this  effort.  MADE-2  was  a  field-scale  natural 
gradient  experiment  performed  in  1990-91  at  Columbus  Air  Force  Base  in  Columbus, 
Mississippi.  A  MADE-2  database  has  been  prepared  by  Boggs  and  others  (1993a)  and 
analyses  have  been  published  by  Boggs  and  others  (1993b)  and  by  Stauffer  and 
others  (1994)  . 

The  MADE-2  test  site  was  an  area  about  300  m  x  200  m  with  about  2  m  of  relief. 
It  was  covered  primarily  by  weeds  and  brush,  and  contained  no  streams  or  ponds. 
The  10  m  to  15  m  thick  upper  layer  of  soil  was  a  shallow  alluvial  terrace 
containing  an  unconfined  aquifer.  This  was  bounded  below  by  an  aquitard  of 
marine  silt  and  clay  (Boggs,  Young,  Benton,  and  Chung;  1990).  The  aquifer  soil 
was  classified  as  poorly  sorted  to  well  sorted  sandy  gravel  and  gravelly  sand 
with  minor  amounts  of  silt  and  clay.  The  aquifer  was  found  to  consist  of 
irregular  lenses  and  layers  having  typical  horizontal  dimensions  on  the  order  of 
8  m  and  typical  vertical  dimensions  on  the  order  of  1  m. 

The  heterogeneity  of  the  MADE-2  site  was  much  greater  than  that  of  other  reported 
natural  gradient  macrodispersion  experiments.  Measurements  using  the  borehole 
flowmeter  method  showed  hydraulic  conductivity  variations  of  up  to  four  orders 
of  magnitude  in  individual  profiles.  Rehfelt,  Boggs,  and  Gelhar  (1992)  found 
that  the  variance  of  the  natural  logarithm  of  the  hydraulic  conductivity  was  at 
least  an  order  of  magnitude  larger  at  Columbus  than  at  Borden,  Twin  Lakes,  or 
Cape  Cod.  The  horizontal  and  vertical  correlation  scales  for  hydraulic 
conductivity  were  also  larger  by  factors  of  1.75  or  more. 

MADE-2  focused  on  the  fate  and  transport  of  dissolved  organic  chemicals  of  the 
types  found  in  jet  fuels  and  solvents.  A  volume  of  9.7  m3  of  tracer  solution  was 
injected  at  a  constant  rate  for  48.5  hours  through  5  wells  spaced  1  m  apart.  The 
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solution  contained  tritiated  water  (an  essentially  passive  tracer) ,  benzene, 
naphthalene,  p-xylene,  and  o-dichlorobenzene .  The  spread  of  the  plume  in  three 
dimensions  was  monitored  for  15  months  by  analyzing  water  samples  drawn  from  up 
to  328  multilevel  sampling  wells  (at  up  to  30  depths  per  well)  and  56  BarCad 
positive  displacement  samplers.  Five  comprehensive  sets  of  water  samples  (called 
snapshots)  were  obtained  at  intervals  of  about  100  days.  Plots  of  concentration 
contours  in  horizontal  planes  showed  that  the  tritium  plume  spread  in  an 
essentially  linear  fashion  with  a  tendency  to  bend  toward  the  northeast.  The 
vertical  structure  along  the  plume  axis  was  complex. 

Boggs  and  others  (1993b) ,  based  on  numerical  integration  of  the  tritium 
concentrations,  found  ratios  of  observed  mass  to  injected  mass  in  the  first  four 
snapshots  of  1.52,  1.05,  0.98,  and  0.77,  respectively.  The  52%  overestimate  in 
the  initial  snapshot  was  attributed  to  preferential  sampling  from  more  permeable 
zones  and  to  vertical  interconnections  between  sampling  points.  The  23% 
underestimate  in  snapshot  4  was  partially  due  to  the  motion  of  the  plume's 
leading  edge  past  the  farthest  downstream  sampling  points.  Snaphot  5  was  not 
intended  to  define  the  entire  plume. 

Our  objective  in  the  1994  Summer  Research  Program  was  to  obtain  improved 
simulations  of  the  MADE -2  tritium  plume  using  public  domain  computer  codes  for 
groundwater  flow  and  contaminant  transport.  The  present  work  is  an  extension  of 
the  senior  author's  previous  efforts  as  an  AFOSR  Summer  Faculty  Fellow  (Gray, 
1992;  1993) . 

FLOW  MODELING 


In  accord  with  most  groundwater  studies,  in  the  present  work  the  effects  of 
density  variations  are  assumed  to  be  negligible,  so  that  the  flow  equation  can 
be  solved  without  knowing  the  concentration  field.  The  resulting  velocity  field 
is  input  to  the  transport  equation,  which  is  then  solved  for  the  concentrations. 
These  calculations  were  performed  using  computer  programs  MODFLOW  for  the  flow 
problem  and  MT3D  for  the  transport  problem.  Many  other  programs  were  used  to 
prepare  input  files  or  to  analyze  results.  Unless  noted  otherwise,  these  were 
written  by  the  authors  of  this  report  in  FORTRAN  77. 

MODFLOW  (McDonald  and  Harbaugh,  1988)  is  a  U.  S.  Geological  Survey  (USGS)  public 
domain  FORTRAN  77  program  for  the  solution  of  the  groundwater  flow  equation.  The 
program's  name  refers  to  its  modular  structure  which  facilitates  the  insertion 
of  new  subroutine  packages  to  handle  specific  tasks.  The  version  used  here, 
MODFLOW/mt,  was  obtained  from  Dr.  Chunmiao  Zheng,  the  author  of  MT3D,  and 


4 


incorporated  several  new  subroutine  packages  which  are  described  below. 
Flexibility,  robustness,  clarity  of  coding,  and  outstanding  documentation  all 
contributed  to  the  selection  of  MODFLOW  for  this  project. 

The  basic  MODFLOW  program  solves  a  block  centered  finite  difference  approximation 
to  the  groundwater  flow  equation  on  a  variable  cell  size,  three  dimensional 
rectangular  grid.  MODFLOW  allows  for  anisotropy  so  long  as  the  grid  axes  are 
aligned  with  the  principal  directions  of  hydraulic  conductivity.  It  can  solve 
either  steady  or  transient  cases  and  provides  options  for  recharge,  wells,  and 
other  hydrologic  features.  Both  confined  and  unconfined  aquifers  can  be  modeled. 
The  original  block  centered  flow  package  (BCF1)  allowed  the  dewatering  of  layers 
during  periods  of  water  table  decline,  but  could  not  handle  rewetting  due  to  a 
rising  water  table.  This  was  an  important  limitation  in  modeling  MADE -2  due  to 
the  pronounced  water  table  fluctuations  which  were  observed.  The  version  used 
here  incorporated  BCF2  (McDonald,  Harbaugh,  Orr,  and  Ackerman;  1991) ,  a  newer 
package  which  allows  rewetting.  The  present  MODFLOW  also  incorporated  PCG2 
(Hill,  1990)  ,  a  preconditioned  conjugate  gradient  solver;  LKMT18 ,  which  generates 
output  files  in  a  format  suitable  for  input  to  MT3D;  and  STR1,  a  stream 
interaction  package  which  was  not  used. 

The  user  of  MODFLOW  must  input  the  grid  geometry,  boundary  and  initial 
conditions,  values  related  to  the  principal  hydraulic  conductivities  for  each 
cell,  storage  coefficients  for  each  cell,  and  source  parameters. 

The  definition  of  a  suitable  computational  grid  is  the  first  step  in  applying 
MODFLOW.  In  view  of  the  heterogeneity  of  the  site  and  the  nature  of  the  plume, 
a  uniform  three  dimensional  grid  was  selected.  As  in  Gray  (1993)  ,  the  grid 
consists  of  9  layers,  each  containing  66  rows  and  21  columns  of  5  m  x  5  m  cells. 
The  horizontal  grid  is  identical  to  that  of  Gray  (1993)  with  the  105  m  and  330 
to  sides  parallel  to  the  x  and  y  axes  of  the  MADE-2  coordinate  system, 
respectively.  The  origin  of  the  MADE-2  coordinate  system  is  at  the  center  of  the 
cell  which  contains  all  5  injection  wells  (row  61,  column  11)  .  In  terms  of  MADE- 
2  coordinates,  the  domain  extends  from  -52.5  m  to  +52.5  m  in  the  x  direction  and 
from  -27.5  m  to  +302.5  m  in  the  y  direction. 

One  of  the  most  critical  steps  in  the  development  of  a  numerical  model  is 
geostatistical  analysis,  the  process  by  which  a  relatively  small  number  of 
irregularly  spaced  observations  of  some  variable  are  used  to  assign  values  at  the 
relatively  large  number  of  regularly  spaced  computational  nodes.  Gray  (1993) 
used  the  commercial  program  SURFER  for  this  task.  In  the  present  study  the 
public  domain  software  package  Geo-EAS  Version  1.2.1  (Englund  and  Sparks,  1991) 
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was  employed.  Geo-EAS  is  a  menu  driven  personal  computer  program  developed  by 
the  Environmental  Protection  Agency  (EPA)  primarily  to  perform  two  dimensional 
kriging.  Geo-EAS  allows  the  user  to  closely  control  most  aspects  of  the  kriging 
process,  including  the  selection  of  linear,  spherical,  exponential,  or  Gaussian 
variograms.  The  program  can  also  calculate  descriptive  statistics  and  produce 
two  dimensional  contour  plots.  In  comparison  with  SURFER  Version  4,  Geo-EAS  is 
less  polished,  has  inferior  graphics,  and  has  more  glitches,  e.g.  the  Gaussian 
variogram  doesn't  always  work.  On  the  other  hand,  Geo-EAS  is  more  flexible  and 
is  much  less  of  a  black  box.  In  this  study  all  kriging  was  done  using  Geo-EAS, 
but  most  of  the  final  contour  plots  were  made  using  SURFER  Version  4. 

Geological  logs  from  32  locations  scattered  over  and  near  the  site  were  analyzed 
to  determine  the  vertical  boundaries  of  the  aquifer.  Program  XLTOGE  was  written 
to  reformat  the  measured  ground  surface  and  aquitard  top  elevations  for  input  to 
Geo-EAS.  These  data  were  kriged  using  a  linear  variogram  for  the  ground  surface 
elevation  and  a  spherical  variogram  for  the  aquifer  bottom  elevation.  The  ground 
surface  elevation  was  estimated  to  vary  from  64.68  m  to  65.99  m,  and  the  aquifer 
bottom  was  estimated  to  range  from  49.90  m  to  55.51  m  MSL. 

The  rewetting  capability  of  the  BCF2  package  allowed  for  a  more  efficient 
vertical  grid  spacing  that  had  been  used  previously.  In  Gray  (1993)  ,  the 
computational  domain  was  bounded  below  by  an  impermeable  plane  at  51.0  m,  and  the 
lower  8  layers  were  each  1  m  thick.  The  top  layer,  with  a  base  at  59.0  m,  had 
an  upper  boundary  which  fluctuated  with  the  water  table.  As  the  observed  water 
table  reached  its  peak  in  May  1991,  cells  in  the  top  layer  were  up  to  6 . 1  m 
thick.  This  was  undesirable  from  the  standpoint  of  accuracy,  but  was  necessary 
because  BCF1  required  the  lower  boundary  of  the  top  layer  to  be  low  enough  to 
guarantee  against  dewatering. 

In  the  present  grid  the  base  of  the  upper  layer  is  at  63.0  m,  so  that  its 
saturated  thickness  should  never  exceed  2.1  m.  The  next  seven  layers  are  each 
1  m  thick.  The  top  of  the  lowest  layer  is  at  56.0  m,  and  its  impermeable  bottom 
varies  to  match  the  top  of  the  aquitard.  The  thickness  of  the  lowest  layer 
ranges  from  0.49  m  to  6.10  m  with  a  mean  of  3.31  m.  In  terms  of  MODFLOW 
classification,  layer  1  is  unconfined,  layers  2  through  7  are  fully  convertible 
(LAYCON  =  3) ,  and  layers  8  and  9  are  confined. 

There  were  82  piezometers  scattered  irregularly  over  and  near  the  computational 
domain.  Heads  were  recorded  continuously  in  16  piezometers.  There  were  also  17 
manual  piezometer  surveys  conducted  at  intervals  of  about  one  month  and  typically 
covering  45  piezometers.  The  continuous  and  survey  observations  showed  good 
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agreement.  From  the  first  observations,  about  1  week  before  injection,  until 
about  180  days  after  injection,  heads  declined  smoothly  less  than  1  m.  After 
that  date  heads  underwent  larger  and  more  erratic  changes.  These  results  showed 
that  a  transient  model  was  essential. 

The  piezometric  heads  from  the  monthly  surveys  were  needed  to  establish  the 
initial  head  at  each  node,  as  well  as  the  head  at  each  boundary  node  as  a 
function  of  time.  Using  SURFER,  Gray  (1993)  kriged  using  all  of  the  available 
heads,  pooling  all  depths  and  including  piezometers  which  were  far  from  the 
computational  domain.  The  results  were  assigned  as  initial  and  boundary 
conditions  to  all  layers,  i.e.  there  was  no  variation  of  head  with  depth.  The 
numerical  solutions  obtained  with  these  conditions  showed  heads  which  dropped 
toward  the  northwest  corner  of  the  grid,  suggesting  that  the  plume  should  bend 
toward  the  northwest.  As  the  observations  showed  the  plume  bending  toward  the 
northeast,  it  was  important  to  be  more  careful  in  translating  the  observed  heads 
into  initial  and  boundary  conditions. 

The  commercial  spreadsheet  Quattro  Pro  for  Windows  was  used  to  examine  the 
distribution  of  the  piezometer  screen  midpoint  elevations.  It  was  noticed  that 
most  were  close  to  either  60.5  m  or  56.0  m.  Geo-EAS  was  used  to  reject 
piezometers  which  were  not  close  to  these  elevations  or  were  too  far  outside  the 
computational  domain.  The  pizometers  selected  for  kriging  consisted  of  an  upper 
set  of  15  whose  screen  midpoints  ranged  from  59.76  m  to  61.22  m  with  a  mean  of 
60.55  m,  and  a  lower  set  of  23  whose  elevations  were  between  55.51  m  and  56.71 
m  with  a  mean  of  55.95  m.  Figure  1  shows  that  the  coverage  of  the  (plan)  north 
end  of  the  computational  domain  was  sparse  at  both  levels. 

MADETOGE  was  written  to  segregate  the  monthly  piezometer  survey  data  into  upper 
and  lower  piezometer  files.  These  files  were  kriged  with  linear  variograms  using 
Geo-EAS.  Figure  2  shows  the  results  for  the  upper  and  lower  piezometer  sets  for 
the  survey  of  March  8,  1991.  In  almost  every  survey  the  heads  at  both  levels 
decline  toward  the  northwest.  The  upper  level  heads  were  assigned  to  layers  1 
through  4,  and  the  lower  level  heads  to  layers  8  and  9.  Heads  were  specified  for 
layers  5,  6,  and  7  by  linear  interpolation.  Program  BASMAKER  wrote  the  MODFLOW 
Basic  package  input  file  which  included  the  initial  heads  at  every  node.  Program 
GHBMAKER  created  the  input  file  for  the  MODFLOW  General  Head  Boundary  package. 
The  function  of  this  package  was  to  maintain  specified  heads  at  every  boundary 
node  (Dirichlet  boundary  conditions) . 

The  net  recharge  was  the  difference  between  precipitation  and 
evapotranspiration.  Daily  precipitation  and  temperature  data  were  measured  at 
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the  CAFB  weather  station,  less  than  2  km  from  the  test  site.  Daily  pan 
evaporation  data  from  State  University,  about  35  km  distant,  were  supplied  by 
State  Climatologist  Dr.  C.  L.  Wax.  Missing  evaporation  data  were  estimated  from 
the  daily  maximum  temperatures  using  the  empirical  equation  of  Pote  and  Wax 
(1986).  Based  on  the  recommendation  of  Dr.  Wax,  a  pan  coefficient  of  0.8  was 
used  to  estimate  the  evapotranspiration . 

The  17  piezometer  surveys  and  the  two  day  injection  period  were  used  to  define 
18  stress  periods  during  which  all  boundary  conditions  and  water  sources  were 
constant.  These  were  the  same  periods  used  by  Gray  (1993).  Except  for  the 
injection  period,  the  stress  periods  were  approximately  centered  on  the  survey 
dates.  The  recharge  rates  were  the  averages  of  the  daily  values.  Table  1 
defines  the  stress  periods  used  in  MODFLOW.  The  injection  occurred  at  a  rate  of 
4.85  m3/ day  on  simulation  days  15  and  16  at  row  61,  column  11,  and  layer  7.  A 
constant  time  step  of  2  days  was  used  in  all  the  MODFLOW  simulations. 


Table  1.  Stress  periods  and  recharge  rates  used  in  MADE * 2  simulations. 


stress 

period 


starting 

date 


starting 
s im .  day 
number 


1 


period 

length 

[days] 


14 


survey 

date 


June  19 


survey  recharge 

sim.  day  rate 


number 


8 


13 

Apr.  22 

14 

May  16 

15 

June  3 

16 

June  2  7 

17 

July  31 

18 

Sept.  1 

last  day 

Sept.  22 

Mar .  8 


Ap: 


[m/day] 


-*0.00313 


-0 . 00478 


“0 . 00148 


-0.00409 


-0.00286 


-0 . 00107 


+0 . 00071 


+0 . 00942 


+0 . 00387 


+0 . 00809 


+0 . 00114 


+0 . 00794 


+0 . 01022 


+0.00357 


+0.00046 


-0 . 00273 


-0 . 00159 


-0 . 00384 
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Vertical  profiles  of  horizontal  hydraulic  conductivity  were  measured  in  67  wells 
scattered  in  and  around  the  computational  domain.  The  data  were  measured  over 
successive  15  cm  layers  using  a  borehole  flowmeter.  The  gaps  where  the  well 
screens  were  jointed  were  filled  in  with  the  values  immediately  above  and  below. 
The  height  profiled  and  the  layer  boundaries  varied  from  well  to  well. 

KAVG94  was  written  to  relate  these  profiles  to  the  grid  layers.  The  tops  of  the 
profiles  varied  from  57.62  m  to  62.68  m.  The  program  extended  each  profile  up 
to  64.0  m  using  the  conductivity  at  the  top  of  the  profile.  The  lowest  points 
varied  from  51.88  m  to  56.22  m.  Profiles  were  extended  down  to  56.0  m  or  the 
next  lower  integer  elevation  using  the  conductivity  at  the  lowest  point.  The 
extended  profiles  were  averaged  arithmetically  over  each  MODFLOW  layer  to 
generate  horizontal  conductivities.  With  the  assumption  that  each  15  cm  slice 
of  material  was  isotropic,  the  extended  profiles  were  averaged  harmonically 
between  the  midpoints  of  the  MODFLOW  layers  to  generate  vertical  leakances. 
Leakance  is  the  vertical  conductivity  divided  by  the  thickness  between  adjacent 
nodes.  Due  to  the  variable  thickness  of  layer  9,  the  leakance  between  layers  8 
and  9  was  calculated  for  the  interval  from  56.5  m  to  55.5  m  rather  than  to  the 
actual  midpoint  of  the  lowest  cells.  Exceptions  occurred  at  wells  K-2,  K-26,  and 
K-28  where  the  profiles  ended  at  56.0  m. 

The  next  task  was  to  interpolate  and  extrapolate  the  averaged  profiles 
horizontally  so  as  to  obtain  the  horizontal  conductivity  and  vertical  leakance 
at  each  node  of  the  computational  grid.  The  averaged  profiles  were  log 
transformed  using  KA2LOG,  kriged  with  Geo-EAS,  and  transformed  back  by  DLOGFILE . 
The  log  transformation  was  necessary  to  avoid  negative  values  in  the  kriging 
process.  Spherical  or  exponential  variograms  were  used.  Program  BCF2MAKER  was 
written  to  format  the  conductivity  values  for  input  to  the  MODFLOW  BCF2  package. 

During  execution,  MODFLOW  calculates  the  transmissivity  of  the  cells  which  are 
partially  saturated  by  multiplying  the  horizontal  conductivity  of  the  cell  by  its 
saturated  depth.  Since  the  horizontal  hydraulic  conductivity  represents  an 
average  over  the  entire  cell  thickness,  this  is  correct  only  if  the  cell  is  truly 
homogeneous.  The  vertical  leakance  is  treated  as  a  constant  as  long  as  a  cell 
contains  water,  even  though  it  represents  an  average  over  the  full  region  between 
nodes.  This  is  not  correct  either. 

Little  was  known  about  the  storage  coefficients.  A  specific  yield  of  0.1  was 
measured  in  a  single  traditional  pump  test  (AT- 2)  (Boggs,  Young,  Benton,  and 
Chung;  1990) .  No  measurements  of  specific  storage  were  made,  so  a  confined 
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storage  coefficient  base  value  of  0.0001  was  assumed,  based  on  textbook  values 
for  specific  storage  in  sand  and  sandy  gravel  (Anderson  and  Woessner,  1992) .  In 
view  of  the  great  uncertainty  of  these  parameters,  simulations  were  run  with 
higher  and  lower  values  in  order  to  investigate  the  sensitivity  of  the  results. 
In  each  simulation,  the  storage  coefficients  were  constant  throughout  the  grid. 
In  reality,  great  variability  is  expected;  but  there  was  no  defensible  way  to 
account  for  this  on  the  basis  of  the  available  data. 


The  468  day  experiment  was  simulated  on  a  Sun  Sparcstation  2  using  the  PCG2 
solver.  In  spite  of  the  rather  severe  vertical  motion  of  the  water  table, 
MODFLOW  performed  reliably.  Table  2  lists  the  differences  among  the  five  cases 
which  were  computed. 


Table  2.  MODFLOW  simulation  summary. 


Case 

RELAX 

WETDRY 

[meters] 

specific 

yield 

confined 
storage 
coef . 

run  time 
[min.  ] 

final 

volume  ! 

error 

1 

0 . 98 

-0.1 

0.1 

0 . 0001 

60 

-0.25% 

2 

H 

O 

O 

-0.1 

0 . 1 

0 . 0001 

unknown 

-0.24% 

3 

CD 

Ch 

o 

-0 . 01 

0 . 1 

0.0001 

72 

-0.25% 

4 

o 

VD 

00 

-0.1 

0.2 

0.0005 

94 

-1.52% 

5 

0 . 98 

-0 . 1 

0.05 

0 . 00005 

58 

-0 . 23% 

Taking  Case  1  as  the  base  case,  Case  2  tests  the  effect  of  increasing  RELAX,  a 
convergence  parameter  in  the  PCG2  solver  package.  This  variation  left  the 
solution  virtually  unchanged.  Case  3  examines  the  effect  of  reducing  WETDRY,  a 
parameter  in  the  BCF2  package  which  controls  cell  rewetting.  The  negative  sign 
indicates  that  the  rewetting  of  cell  x  depends  on  the  head  in  the  cell  below. 
The  absolute  value  of  WETDRY  is  the  amount  by  which  the  head  in  the  cell  below 
must  exceed  the  bottom  elevation  of  cell  x  before  it  rewets.  Case  3  results  were 
virtually  identical  with  Case  1.  A  positive  value  of  WETDRY  makes  rewetting 
depend  on  the  heads  in  the  four  horizontally  adjacent  cells.  Runs  with  positive 
values  of  WETDRY  invariably  failed  to  converge. 

Cases  4  and  5  varied  the  storage  coefficient  values.  It  can  be  seen  that 
increasing  the  storage  coefficients  increases  the  volumetric  discrepancy  and  the 
run  time.  The  effects  on  the  nature  of  the  solution  are  discussed  further  below, 
but  they  have  not  yet  been  fully  assessed. 

Figure  3  presents  the  Case  1  head  contours  on  simulation  day  270  (March  8,  1991) 
in  layers  4  and  9.  Compared  with  the  kriged  distributions  for  the  upper  and 
lower  piezometers  on  the  same  day  shown  in  Figure  2,  it  can  be  seen  that  the  head 
distributions  are  both  qualitatively  and  quantitatively  similar.  In  both  the 
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predicted  and  observed  cases,  the  flow  is  downward.  The  tendency  for  the  heads 
to  decline  toward  the  northwest  is  evident  in  this  figure  and  throughout  the 
simulation. 

In  order  to  obtain  a  numerical  measure  of  agreement,  the  simulated  heads  were 
compared  to  the  continuous  head  observations.  Program  WELLGRPH  was  written  to 
extract  from  the  MODFLOW  binary  output  file  the  head  time  series  for  those  cells 
which  contained  continuously  monitored  piezometers.  The  continuous  piezometer 
records  show  erratic  day  to  day  variations  which  cannot  be  predicted  by  a  model 
whose  boundary  conditions  change  only  16  times  in  468  days.  To  provide  a 
reasonable  basis  of  comparison,  the  daily  observed  heads  were  averaged  over  each 
stress  period  by  program  HYDROGRA.  Figure  4  compares  the  Case  1  predictions  to 
the  observed  (averaged)  heads  at  two  piezometers  with  the  same  horizontal 
position.  The  simulated  results  adjust  rapidly  to  the  boundary  conditions  for 
each  stress  period.  The  model  results  are  better  at  the  upper  level  (P55a)  than 
at  the  lower  level  (P55b) ,  where  the  model  overpredicts  markedly  in  stress 
periods  9,  11,  and  13. 

The  averaged  observations  were  subtracted  from  the  unaveraged  MODFLOW  heads  and 
the  maximum,  minimum,  and  root  mean  square  (rms)  differences  were  summarized  in 
Table  3.  Case  5,  with  the  smallest  storage  coefficients,  gives  the  best  overall 
accuracy.  Case  4  has  the  greatest  excursions  from  the  observations,  yet  its  rms 
deviation  is  smaller  than  Case  1.  Although  the  ability  of  the  model  to  reproduce 
the  observations  is  imperfect,  it  is  hard  to  see  how  the  model  could  be  improved 
given  the  limitations  of  the  data  base. 


Table  3.  Deviation  of  MODFLOW  heads  from  continuous  observations  [meters]. 


max. 

max. 

max. 

rms 

rms 

rms 

Well 

Case  1 

Case  4 

Case  5 

Case  1 

Case  4 

Case  5 

Case  1 

Case  4 

Case  5 

P53a 

-0.65 

-1.32 

-0.57 

0.74 

0.51 

0.14 

0.329  ; 

0.228 

0.194 

P54a 

gym 

-0.37 

0.39 

0.58 

0.30 

0.143 

0.165 

0.136 

P54b 

-0.42 

-0.78 

-0.17 

0.43 

0.52 

0.43 

0.147 

0.159 

0.143 

P55a 

-0.53 

-0.80 

-0.37 

0.44 

0.50 

0.44 

0.199 

0.204 

0.199 

P55b 

mm 

-0.44 

+0.01 

1.01 

1.01 

1.01 

0.374 

0.374 

0.374 

P60a 

-0.51 

-1.51 

0.30 

0.38 

0.30 

0.188 

0.188 

0.188 

P61a 

I 

-0.40 

gum 

-0.40 

0.36 

0.36 

0.36 

0.188 

0.188 

0.188 

P61b 

-0.39 

0.23 

0.23 

0.23 

0.154 

0.154 

0.154 

average 

-0.69 

-0.35 

0.49 

0.51 

0.40 

0.215 

0.208 

0.197 
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TRANSPORT  MODELING 


MT3D  is  a  public  domain  program  developed  for  the  EPA  to  solve  the  three 
dimensional  groundwater  transport  equation  for  dissolved  contaminants  (Zheng, 
1990)  .  MT3D  is  coded  in  Fortran  77  and  uses  the  same  modular  structure  as 
MODFLOW.  In  fact,  MT3D  accepts  as  input  the  head  and  flux  distributions  computed 
by  MODFLOW  (or  similar  flow  models)  .  MT3D  then  predicts  the  concentration  field 
of  a  single  contaminant  which  undergoes  advection,  dispersion,  and  chemical 
reactions .  The  program  provides  for  various  types  of  point  and  area  sources  and 
sinks  including  wells,  recharge,  and  flows  through  the  domain  boundaries.  MT3D 
Version  1.80  was  used  in  this  study. 

Because  of  the  computational  difficulties  of  numerical  dispersion  and  oscillation 
in  advection-dominated  flows,  MT3D  incorporates  four  options  for  calculating  the 
advection  term.  The  Method  of  Characteristics  (MOC)  tracks  a  large  number  of 
imaginary  tracer  particles  forward  in  time.  The  Modified  Method  of 
Characteristics  (MMOC)  tracks  particles  located  at  the  cell  nodes  backward  in 
time.  The  MMOC  requires  much  less  computation  than  the  MOC,  but  it  is  not  as 
effective  in  eliminating  artificial  dispersion,  especially  near  sharp  fronts. 
The  Hybrid  Method  of  Characteristics  (HMOC)  uses  the  MOC  near  sharp  concentration 
gradients  and  the  MMOC  in  the  remainder  of  the  domain.  An  Eulerian  Upstream 
Differencing  (UD)  option  is  provided  for  problems  in  which  advection  does  not 
dominate . 

The  dispersion  terms  are  computed  using  a  fully  explicit  Eulerian  central 
difference  method.  For  isotropic  media,  the  dispersion  coefficients  are  based 
on  longitudinal  and  transverse  dispersivities .  For  more  complex  situations, 
there  is  an  option  which  distinguishes  horizontal  and  vertical  transverse 
dispersivities.  The  explicit  formulation  reduces  the  memory  needed,  but  requires 
limits  on  the  time  step  to  assure  numerical  stability.  Consequently  each  flow 
model  time  step  may  be  automatically  subdivided  into  several  transport  steps  in 
order  to  maintain  numerical  stability  in  MT3D . 

MT3D  allows  both  equilibrium  sorption  and  first  order  irreversible  rate 
reactions .  Equilibrium  sorption  reactions  transfer  contaminant  between  the 
dissolved  phase  and  the  solid  phase  (which  is  sorbed  to  the  soil  matrix)  at  time 
scales  much  shorter  than  those  of  the  flow.  These  reactions  may  be  described  by 
linear  isotherms  or  nonlinear  isotherms  of  the  Freundlich  or  Langmuir  types.  In 
first  order  irreversible  rate  reactions  the  rate  of  mass  loss  is  linearly 
proportional  to  the  mass  present.  This  class  includes  radioactive  decay  and 
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certain  types  of  biodegradation. 


MT3D  requires  information  beyond  that  needed  for  and  calculated  by  MODFLOW.  A 
porosity  is  needed  for  each  cell  in  order  to  calculate  seepage  velocities,  yet 
porosities  were  measured  in  only  four  core  holes.  The  84  samples  had  a  mean 
porosity  of  0.32,  and  this  value  was  used  for  every  cell.  Based  on  the  MADE -2 
observations  and  an  assumed  two  dimensional  analytical  model  for  the  plume,  Boggs 
and  others  (1993b)  estimated  the  longitudinal  dispersivity  to  be  10  m  and  the 
transverse  horizontal  dispersivity  to  be  less  than  2.2  m.  The  base  values  of 
dispersivity  used  were  10  m  in  the  longitudinal  direction,  1  m  in  the  horizontal 
transverse  direction,  and  0.1  m  in  the  vertical  transverse  direction.  For  the 
purpose  of  calculating  concentrations,  every  wetted  layer  was  assumed  to  have  a 
uniform  thickness  of  1  m,  although  the  actual  thickness  varied  for  the  top  and 
bottom  layers. 

MT3D  was  applied  only  to  the  tritium  plume.  The  molecular  diffusion  coefficient 
of  tritium  in  water,  calculated  using  the  Wilke-Chang  method,  was  multiplied  by 
an  assumed  tortuosity  of  0.25  to  yield  the  value  of  2.16  x  10'4  m2/day  for  the 
molecular  diffusion  coefficient  of  tritium  in  a  saturated  porous  medium.  The 
injected  fluid  had  a  tritium  concentration  of  0.0555  Ci/m3;  and  the  natural 
background,  including  recharge  and  boundary  inflows,  was  set  to  zero.  Water 
leaving  the  domain  carried  the  concentration  of  the  cell  it  last  occupied. 
Sorption  does  not  affect  tritiated  water,  but  tritium  decays  with  a  12.26  year 
half-life . 

The  transport  simulations  attempted,  all  based  on  the  Case  1  MODFLOW  head 
solution,  are  summarized  in  Table  4.  None  are  remotely  satisfactory.  No  run 
extended  beyond  simulation  day  141  because  by  that  time  each  had  experienced  a 
numerical  failure  or  had  been  terminated  because  the  solution  was  unreasonable. 
In  general,  the  run  times  were  inconveniently  long.  The  mass  discrepancies 
appear  either  unacceptably  large  (MOC,  MMOC,  and  HMOC)  or  remarkably  tiny  (UD)  , 
but  the  meaning  of  this  parameter  is  not  clear.  Runs  7  (HMOC)  and  8  (UD) 
predicted  nearly  identical  plumes  even  though  the  mass  discrepancies  were  very 
different . 

Run  3  produced  a  widely  spread  plume  even  though  the  dispersion  package  was 
turned  off.  This  appears  to  be  a  numerical  shortcoming  of  the  MMOC  method 
because  no-dispersion  runs  5  (MOC)  and  6  (HMOC)  predicted  unrealistically  small 
spreads.  All  of  the  no-dispersion  runs  were  free  from  negative  concentrations. 
Run  11  was  a  repetition  of  Run  9  using  double  precision  arithmetic;  the  results 
were  identical.  In  Runs  9  and  11  the  dispersivities  in  the  longitudinal, 


13 


transverse  horizontal,  and  transverse  vertical  directions  were  4.0m,  0.4  m,  and 
0.4  m,  respectively.  Runs  12  (UD)  and  14  (HMOC)  used  dispersivities  in  the 
longitudinal,  transverse  horizontal,  and  transverse  vertical  directions  of  1.0 
m,  0.1m,  and  0.1m,  respectively.  In  Run  13  (UD)  the  dispersivities  were  all  0, 
but  molecular  diffusion  was  active.  In  general,  the  smaller  the  dispersivities, 
the  more  realistic  the  plume  appeared. 


Table  4.  Summary  of  MT3D  simulations. 
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CONCLUSIONS 


1.  Geo-EAS  Version  1.2.1  is  technically  superior  to  SURFER  Version  4.  It 
provides  a  satisfactory  tool  for  exploratory  data  analysis  and  two  dimensional 
kriging.  SURFER  has  better  graphic  capabilities. 

2 .  Three  dimensional  groundwater  flow  simulations  using  MODFLOW  are  practical  and 
consistent.  The  rewetting  capability  of  the  BCF2  package  improves  the  accuracy 
of  simulations  in  which  the  water  table  fluctuates  as  much  as  in  MADE -2 . 

3.  Although  the  flow  model  has  not  been  subjected  to  grid  refinement  or  extensive 
parametric  variation  studies,  the  comparison  between  the  simulated  and  observed 
heads  is  satisfactory.  Given  the  existing  data,  there  is  little  prospect  for 
significant  improvement. 

4 .  The  simulated  head  distribution  suggests  that  the  plume  should  bend  toward 
the  northwest.  The  observations  show  the  plume  bending  toward  the  northeast. 
This  discrepancy  is  probably  due  to  inaccurate  head  boundary  conditions  caused 
by  a  lack  of  piezometers  near  the  northern  end  of  the  grid. 

5 .  We  were  unsuccessful  in  our  attempts  to  simulate  the  spread  of  the  tritium 
plume  using  MT3D .  Further  efforts  to  achieve  complete,  accurate  simulations  of 
the  tritium  plume  should  be  made . 
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Figure  1 . 
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Figure  3.  MODFLOW  Case  1  simulated  heads  for  layers  4  (left)  and  9  (right)  for 
simulation  day  270  (March  8,  1991) .  Heads  are  in  meters. 


